# -*- coding: utf-8 -*-
"""
Created on Fri Oct 29 21:09:51 2021

@author: Decay
"""

import numpy as np

#n=1
def NC_1(a,b,func):
    return 1/2*(b-a)*(func(a)+func(b))
#n=2
def NC_2(a,b,func):
    return 1/6*(b-a)*(func(a)+4*func(a+b)+func(b))
#n=4
def NC_4(a,b,func):
    n=4
    h=(b-a)/4
    f=[]
    for i in range(5):
        f.append(func(a+i*h))
    return 1/90*(7*f[0]+32*f[1]+12*f[2]+32*f[3]+7*f[4])

def func(x):
    if x==0:
        x=1e-20
    return np.sin(x)/x
a=0
b=10

nc1=NC_1(a, b, func)
nc2=NC_2(a, b, func)
nc4=NC_4(a, b, func)

from scipy import integrate
def f(x):
    return np.sin(x)/x
print(integrate.quad(f,a,b))  # quad方法会返回精确的值和误差

print("NC_1 err: %s" %(abs(nc1-integrate.quad(f,a,b)[0])))
print("NC_2 err: %s" %(abs(nc2-integrate.quad(f,a,b)[0])))
print("NC_3 err: %s" %(abs(nc4-integrate.quad(f,a,b)[0])))
